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Abstract 

The paper deals with the interaction between buckhng and resonance instabihties 
of mechanical systems. Taking into account the effect of geometric nonlinearity 
in the equations of motion through the geometric stiffness matrix, the problem is 
reduced to a generalized eigenproblem where both the loading multiplier and the 
natural frequency of the system are unknown. According to this approach, all the 
forms of instabilities intermediate between those of pure buckling and pure forced 
resonance can be investigated. Numerous examples are analyzed, including: discrete 
mechanical systems with one to n degrees of freedom, continuous mechanical systems 
such as oscillating deflected beams subjected to a compressive axial load, as well as 
oscillating beams subjected to lateral-torsional buckling. A general finite element 
procedure is also outlined, with the possibility to apply the proposed approach to 
any general bi- or tri-dimensional framed structure. The proposed results provide a 
new insight in the interpretation of coupled phenomena such as flutter instability 
of long-span or high-rise structures. 

Key words: Buckling; Resonance; Flutter; Discrete systems; Continuous systems; 
Finite elements. 



1 Introduction 



Buckling, resonance and flutter are the main forms of instability of the elastic 
equilibrium of structural systems. In buckling instability, by removing the hy- 
pothesis of small displacements so that the deformed structural configuration 
can be distinguished from the undeformed one, it is possible to show that the 
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solution of an elastic problem can represent a condition of stable, neutral or 
unstable equilibrium, depending on the magnitude of the applied load. As is 
well-known, buckling is usually observed in slender structural elements sub- 
jected to a compressive stress field, such as columns of buildings, machine 
shafts, struts of trusses, thin arches and shells. In some cases, elastic instabil- 
ity can also take place for special loading and geometrical conditions, as for 
example in the lateral torsional buckling of slender beams [1] . 

The phenomenon of resonance is also particularly important in structural en- 
gineering. It represents a form of dynamic instability, which occurs when an 
external periodic frequency matches one of the natural frequencies of vibration 
of the mechanical system. In this case, therefore, structural design deals with 
the determination of such dangerous natural frequencies according to modal 
analysis [2-4]. 

Finally, the phenomenon of flutter is a form of acroelastic instability observed 
in long-span or high-rise structures subjected to wind loads, such as towers, 
tall buildings [5] and suspended [6-8] or cable-stayed bridges [9]. In this case, 
the instability is attributed to motion-induced or self-excited forces, which are 
loads induced or infiuenced by the deformation of the structure itself [10, 11]. 
Forces originated in this way modify the initial deformation of the mechan- 
ical system which, consequently, leads to modified forces, and so on. This 
feed-back mechanism can give rise to an amplification effect on the initial 
deformation, leading to premature failure of the structure. The well-known 
dramatic Tacoma Narrows bridge disaster of 1940 is a famous example of this 
catastrophic interaction, and it is still very much in the public eye today. In 
this field, which is not at present completely understood, aeroelastic instability 
is often considered as the result of the interaction between buckling (static) 
and resonance (dynamic) instabilities. However, only a few theoretical formu- 
lations have been proposed for modelling aerodynamic forces and, in most 
investigations, empirical models are set up in which the parameters related to 
the fluid-structure interaction are established by experiments [12]. 

In the present contribution, we deal with the phenomenon of interaction be- 
tween buckling and resonance instabilities. A state-of-the-art survey of the 
existing Literature shows that this problem has been mainly addressed in the 
field of multi-parameter stability theory [13-20], where the conditions for sta- 
bility of a mechanical system are studied with reference to a perturbation of 
the problem parameters. In this framework, instability domains of some con- 
tinuous oscillatory systems subjected to buckhng loads were provided [17, 18]. 
However, to make the problem analytically treatable, the analysis was mainly 
limited to certain continuous mechanical systems such as defiected beams and 
beams experiencing lateral torsional buckling. More importantly, the concept 
of geometric stiffness matrix, typical of structural engineering approaches, was 
totally neglected. 
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On the other hand, in the field of bridge engineering, the infiuence of the geo- 
metric nonhnearity is usually taken into account by including the contribution 
of the self-induced aerodynamic forces to the applied loads (see e.g. [9,21,22] 
and [23] for a detailed overview of the mathematical methods) . Such noncon- 
servative aerodynamic forces are put in relationship with the displacements 
and the velocities of the points of the mechanical system according to the so- 
called flutter derivatives that have to be experimentally determined in the lab- 
oratories [24]. In these approaches, the onset of flutter instability corresponds 
to the condition of vanishing structural dumping. Therefore, structural damp- 
ing seems to play a fundamental role, although it is practically impossible to 
be analytically evaluated but only experimentally estimated [21]. 



Prom the mathematical point of view, it is important to note that pure buck- 
ling, pure resonance, and also flutter instabilities are usually mathematically 
treated as eigenvalue problems. In this paper, we propose a mathematical the- 
ory for the analysis of buckling and resonance interaction, with the possibility 
to give an insight onto the mechanisms leading to flutter instability. In the 
mathematical treatment, a special focus will be given to the role played by 
the geometric stiffness matrix, which contributes to the reduction of the global 
elastic stiffness due to the effect of the geometric nonhnearity. This represents 
a novelty of our approach with respect to the models available in the Litera- 
ture. As it will be shown in the sequel, the use of the geometric stiffness matrix 
may provide the proper link between multi-parameter stability theory, typical 
of rational mechanics, and the bridge engineering approach, typical of bridge 
engineers. 



According to this formulation, we will demonstrate that the interaction be- 
tween buckling and resonance leads to a generalized eigenvalue problem where 
both the buckling loads and the natural frequencies of the system are unknown 
and represent the eigenvalues. This approach will permit to inspect all the 
forms of structural instability intermediate between pure buckling and pure 
resonance. These limit cases are instead observed either when the dynamics 
of the system is neglected, or when the external buckling forces are equal to 
zero. 



The effectiveness of the proposed methodology will be demonstrated with re- 
spect to not only discrete mechanical systems with one to n degrees of freedom, 
but also for continuous mechanical systems such as oscillating deflected beams 
and beams showing lateral-torsional coupled deformations. Finally, a general 
procedure is established in the framework of the finite element method. 
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Figure 1. Scheme of the first one-degree of freedom system analyzed. 
2 Discrete mechanical systems 

2.1 Discrete mechanical systems with one degree of freedom 



Let us consider the mechanical system shown in Fig.l, consisting of two rigid 
rods connected by an elastic hinge of rotational rigidity k and constrained at 
one end by a hinge and at the other by a roller support. A mass m is placed in 
correspondence of the intermediate elastic hinge and the system is loaded by a 
horizontal axial force A^. Considering the absolute rotation (p of the two arms 
as the generalized coordinate, the total potential energy, W, and the kinetic 
energy, T, of the whole system are: 



k{2^f - 2iV/(l - cosv?) 

2 



-m 



dt 



(/sin 9?) 



1 

+ -m 



dt 



{I — I cos (fi) 



-ml'^^fp'. 



(1) 



The equation of motion can be determined by writing the Lagrange's equation: 



d (dT\ dT 



dW 



(2) 



dt \d(p J d(p d(p 
In the present case, this yields: 

mf(fi = -Akip + 2NI sin (/?, (3) 
which can be suitably linearized in correspondence of = 0: 

ml'^(p = -^kip + 2Nlip. (4) 



Looking for the solution to Eq.(4) in the general form 93 = ^p^e^'^^, where uj 
denotes the natural angular frequency of the system, we obtain the following 
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Figure 2. Scheme of the second one-degree of freedom system analyzed, 
equation which provides the conditions of equihbrium of the system: 

(ik - 2m - u'^mf) ipo = 0. (5) 

A nontrivial solution to Eq.(5) exists if and only if the term in brackets is 
equal to zero. This critical condition corresponding to the bifurcation of the 
equilibrium establishes a one-to-one relationship between the applied axial 
force, N, and the angular frequency, lu: 

2k ml 2 /^x 
JV=-j--^c.^. (6) 

Moreover, Eq.(6) admits two important limit conditions for, respectively, = 
and m = 0. In the former case, Eq.(6) gives the natural angular frequency 



of the system according to pure modal analysis, that is Ui = yAk/mP. In the 
latter, the pure critical Eulerian load is obtained, that is Ni = 2k /I. Dividing 
Eq.(6) by iVi, we obtain the following relationship between and a; in a 
nondimensional form: 



As a second example, let us consider the mechanical system shown in Fig. 2, 
consisting in two rigid rods on three supports, of which the intermediate one is 
assumed to be elastically compliant with rigidity k. As in the previous case, a 
mass m is placed in correspondence of the intermediate hinge and the system 
is loaded by a horizontal axial force N. Considering the absolute rotation cp 
of the two arms as the generalized coordinate, the total potential energy, W, 
and the kinetic energy, T, of the whole system are: 

W{^) = -k{l sin^f -2NI{1 - cos^), 

\ (8) 

Following the procedure discussed above, we determine the equation of motion 
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by employing the Lagrange's equation (2): 



mP(f — —I sin if{kl cos if — 2N) , 



(9) 



which can be suitably linearized in correspondence of <^ = 0: 



(10) 



Looking for the solution to Eq.(lO) in the general form ip — ipoe^^^, where uj 
denotes the natural angular frequency of the system, we obtain the following 
condition of equilibrium of the system: 



As in the previous example, by setting equal to zero the term in brackets, we 
obtain a one-to-one relationship between the applied axial force, N, and the 
angular frequency, cu: 



This equation admits two important limit conditions for, respectively, N — 
and m = 0. In the former case, Eq.(12) gives the natural angular frequency 



of the system according to pure modal analysis, that is Ui = Jk/m. In the 



latter, the pure critical Eulerian load for buckling instability is obtained, that 
is A^i = kl/2. Dividing Eq.(12) by A^i, we obtain the same relationship between 
the nondimensional terms N/Ni and {uj/ujif as in the previous example (see 



A graphical representation of the condition (7) in Fig.3 shows that the reso- 
nance frequency is a decreasing function of the compressive axial load. This 
demonstrates, for the mechanical systems with a single degree of freedom, 
that the condition of bifurcation of the equilibrium can be reached for a com- 
pressive axial force, A^, lower than the Eulerian buckling load, A^i, provided 
that the system is subjected to an external excitation with frequency uj given 
by Eq.(7). Conversely, failure due to resonance can take place for a; < a;i, 
provided that the system is loaded by an axial force N given by Eq.(7). 

Finally, the issue of stability or instability of the mechanical system in the 
correspondence of the bifurcation point can be discussed as in the static case, 
i.e. by evaluating the higher order derivatives of the total potential energy W . 




(11) 




Eq.(7)). 
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Figure 3. Nondimensional frequency vs. nondimensional axial force for the single 
degree of freedom systems. 

2.2 Discrete mechanical systems with n degrees of freedom 

Let us consider the mechanical system with two degrees of freedom shown 
in Fig.4, consisting of three rigid rods connected by two elastic hinges of 
rotational rigidity /c, and constrained at one end by a hinge and at the other 
by a roller support. A mass m is placed in correspondence of the intermediate 
elastic hinges and the system is loaded by a horizontal axial force N . Assuming 
the vertical displacements x-i and of the elastic hinges as the generahzed 
coordinates, the total potential energy, W, and the kinetic energy, T, of the 
whole system are given by: 



Performing a Taylor series expansion of Eq.(13) about the origin, and assuming 
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Figure 4. Scheme of the two-degrees of freedom system analyzed. 
Xi/l < 1/10 and X2/I < 1/10, we obtain: 

W(xi, X2) = ^ (5xj + 5x1 - &X1X2) - y (2;? + ^2 - X1X2) , 
T(a;i,a;2) = ^mxi^ + ]^mx2^. 



(14) 



The equations of motion are identified by considering the Lagrange's equa- 
tions: 



dt \ dxj 



dT 
dxi 



dW 
dxi ' 



1,2. 



(15) 



In matrix form, they are: 



m 


1 




w 


m 









4k 

Ak 5k 
1? ¥ 



Xi 
X2 



1 2' 

7 7 



X2 



(16) 



Looking for the solution to Eq.(16) in the general form {g} = {qo}e^^^^ where 
uj denotes the natural angular frequency of the system, we obtain the following 
equation, written in symbolic form: 

(-u;2[M] + [K] - N[K,]) {go} = {0}, (17) 

where [M], [K] and [Kq] denote, respectively, the mass matrix, the elastic 
stiffness matrix and the geometric stiffness matrix of the mechanical system. 
Their expressions can be simply obtained by comparying Eq.(17) with Eq.(16). 

A nontrivial solution to Eq.(17) exists if and only if the determinant of the re- 
sultant coefficient matrix of the vector {go} vanishes. This yields the following 
generalized eigenvalue problem: 

det [[K] - N[Kg\ - cu2[M]) = 0, (18) 
where N and o;^ represent the eigenvalues. For this example, Eq.(18) provides 
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the following relationships between the eigenvalues u and N: 



2_ k N 
mP ml ' 

As limit cases, if m = 0, then we obtain the Eulerian buckling loads: 

N, = J, (20a) 
= y , (20b) 
whereas, if N — 0, we obtain the natural frequencies of the system: 



As far as the eigenvectors are concerned, the system (17) yields the eigenvec- 
tors corresponding, respectively, to the eigenfrequencies (19a) and (19b) as 
functions of N: 

^k/l - N 

Xi = —TT, a;2, 22a 



Dividing Eqs.(19a) and (19b) by uj\^ we derive the following nondimensional 
relationships between the eigenvalues: 

-) =l-(-)> (23a) 

(23b) 



- 



Ni \Ni 

In analogy with the results for the single degree of freedom systems, a graphi- 
cal representation of Eqs.(23a) and (23b) is provided in Fig. 5. We notice that 
both the eigenfrequencies are decreasing functions of the compressive axial 
load. Starting from N = 0, bifurcation of the equilibrium would correspond 
to pure resonance instability. Entering the diagram with a value of the nondi- 
mensional compressive axial force in the range < N/Ni < 1, the coordinates 
of the points of the two curves provide the critical eigenfrequencies of the 
mechanical system leading to bifurcation. Axial forces higher than A^i in the 
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Figure 5. Nondimensional frequencies vs. nondimensional axial forces for the two-de- 
grees of freedom system in Fig.4. 



range 1 < N/Ni < N2/N1 can only be experienced if an additional constraint 
is introduced into the system. Moreover, we observe that the apphcd com- 
pressive load influences all the eigenfrequencies. In particular, for the present 
example, the influence of the axial load is greater on the highest frequency 
than on the lower one. 



As a second example of a system with two degrees of freedom, let us examine 
that of Fig.6, which consists of three rigid rods on four supports, of which the 
central ones are assumed to be elastically compliant with rigidity k. A mass m 
is placed in correspondence of the intermediate hinges and the system is loaded 
by a horizontal axial force N. Assuming the vertical displacements Xi and X2 
of the elastic hinges as the generalized coordinates, the total potential energy, 
W, and the kinetic energy, T, of the whole system are given by {xi/l < 1/10 
and X2/I < 1/10): 



W{xi, X2) =^/c (xl + xl^ — Nl 3 — cos ^arcsin 



f . X2\ . 

COS I arcsm ~^ ) ~ cos I arcsm 

N 



X2 — Xl 



■h (x\ + 



2" v~i • -^2) ^ (^1 + 
T{xi,X2) ^-mxi^ + -rnxx" . 



2Cn 



X1X2 



(24) 
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Figure 6. Scheme of the second two-degrees of freedom system analyzed. 
In this case, the Lagrange's equations (15) yield the following matrix form: 











A- 
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> -N 
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Xi 


hi 




m 








k 

















(25) 



Looking for the solution to Eq.(25) in the general form {q\ = {q'oje^'^*, where 
cu denotes the natural angular frequency of the system, we obtain the following 
equation, written in symbolic form: 



(-a;2[M] + [K] - N[K,]) {q^} = {0}, 



(26) 



where [M], [K] and [Kg] denote, respectively, the mass matrix, the elastic 
stiffness matrix and the geometric stiffness matrix of the mechanical system. 
As it can be readily seen, the geometric stiffness matrix for this problem is 
the same as that of the previous example. 

A nontrivial solution to Eq.(26) exists if and only if the determinant of the 
resultant coefficient matrix of the vector {qq} is equal to zero. This yields the 
following generalized eigenvalue problem: 



det ([K] - N[Kg] - u;^[M]) = 0, 



(27) 



where N and uj^ are the eigenvalues of the system. For this example, Eq.(27) 
provides the following relationships between the eigenvalues: 



2 k N 

u = 3—, 

m ml 

2 k N 



m ml 

As limit cases, if m = 0, we obtain the Eulerian buckling loads: 



= -kl, 



N2 = kl, 



(28a) 
(28b) 



(29a) 
(29b) 



11 



1 N2/Ni=3 

Figure 7. Nondimensional frequencies vs. nondimensional axial forces for the two-de- 
grees of freedom system in Fig. 6. 

whereas, if = 0, then we obtain the natural frequencies of the system: 



a;i = a;2 = \ —. (30) 
V m 



As far as the eigenvectors are concerned, the system (26) yields the eigenvec- 
tors corresponding, respectively, to the eigenfrequencies (28a) and (28b), as 
functions of the axial force, A^: 

Nil , , , 



Dividing Eqs.(28a) and (28b) by uj\, we obtain the following nondimensional 
relationships between the eigenvalues: 

5^)^ — ©' 



A graphical representation of Eqs.(32a) and (32b) is provided in Fig. 7. Also 
in this case, both the frequencies are decreasing functions of the compressive 
axial load. However, for the present example, the influence of the axial load is 
greater on the lower frequency of the system than on the higher one. 
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Figure 8. Undeformed and deformed configurations of a deflected beam under com- 
pressive axial force. 

3 Continuous mechanical systems 



3. 1 Oscillations of deflected beams under compressive axial loads 



Let us consider a slender elastic beam of constant cross-section, inextensible 
and not deformable in shear, though deformable in bending, constrained at 
one end by a hinge and at the other by a roller support, loaded by an axial 
force, (see Fig. 8). In this case, with the purpose of analyzing the free flex- 
ural oscillations of the beam, the differential equation of the elastic line with 
second-order effects can be written by replacing the distributed load with the 
force of inertia: 

where El denotes the flexural rigidity of the beam and is its linear density 
(mass per unit length). Equation (33) can be rewritten in the following form: 

where we have set = N/EI. 

Equation (34) is an equation with separable variables, the solution being rep- 
resented as the product of two different functions, each one depending on a 
single variable: 

v{z,t)=7]{z)f{t). (35) 
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Introducing Eq.(35) into Eq.(34), leads: 



d7^+^d;^^+^^d^^°- ^''^ 

Dividing Eq.(36) by the product rjf, we find: 

dV A , 

dt2 _ d^4 ^ dz^ _ , ,2 

where u"^ represents a positive constant, the left and the right hand-sides 
of Eq.(37) being at the most functions of the time t and the coordinate z, 
respectively. From Eq.(37) there follow two ordinary differential equations: 



dt2 



+ a;V = 0, (38a) 



with 

(-) 

Whereas Eq.(38a) is the equation of the harmonic oscillator, with the well- 
known complete integral 

f{t) ^ Acosut + Bsinut, (40) 
Eq.(38b) has the following complete integral 

r]{z) = Ce^'' + De^^' + Ee-^'' + Fe-^^', (41) 
where Ai and A2 are functions of a and /3: 



A,, . ^-J1±^E±^. (42) 

As in the modal analysis, the constants A and B can be determined on the 
basis of the initial conditions, while the constants C, D, E and F can be 
determined by imposing the boundary conditions. As it will be shown in the 
sequel, for a given value of /3, the parameters uj and a can be determined by 
solving a generalized eigenvalue problem resulting from the imposition of the 
boundary conditions. From the mathematical point of view, this eigenvalue 
problem is analogous to that shown for the discrete systems. On the other 
hand, since we are considering a continuous mechanical system having infinite 
degrees of freedom, we shall obtain an infinite number of eigenvalues uji and 
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ai, just as also an infinite number of eigenfunctions fi and rji. Tlie complete 
integral of the differential equation (33) may therefore be given the following 
form, according to the Principle of Superposition: 



v{z,t) = J2vi{z)Mt), (43) 

i=l 

with: 

fi{t) = Ai cos uJit + Bi sin ujit, (44a) 
ri{z) = Qe^'^' + Die^^^^ + Eie'^^'^ + Fie'^^'^ (44b) 

It is important to remark that the eigenfunctions rji are still orthonormal func- 
tions, as in the classical modal analysis (sec the mathematical demonstration 
reported in the Appendix). This permits to determine the constants and 
Bi in Eq.(44a) via the initial conditions (see also [1], Pag. 315): 

v(z^O) ^ Vo{z), (45a) 
dv 

-{z = 0)=Mz). (45b) 



As regards the boundary conditions, let us consider as an example a beam 
supported at both ends, of length /: 



'r^(O) 


= 0, 


1111 




c 







^ V'(o) 


= 0, 


^1 -^2 -^1 -^2 




D 













< 




> = < 
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' riil) 


= 0, 


Q^il c^'-^^ c c 




E 









= 0, 


X2^x,i xle^2i xje-^ii X^q-^^i 




F 








For a nontrivial solution to the system in Eq.(46), the determinant of the co- 
efficient matrix has to vanish. The resulting eigenequation permits, for each 
given value of the parameter P, to determine the eigenvalues ai of the sys- 
tem. Finally, the corresponding natural eigenfrequencies cji can be obtained 
by inverting Eq.(39). 

As an illustrative example, the first three nondimensional frequencies of the 
simply supported beam shown in Fig. 8 are reported in Fig. 9 as functions of the 
applied nondimensional axial force. Parameters cui and Ni denote, respectively, 
the i-th frequency of the system determined according to modal analysis and 
the i-th buckling load determined according to the Euler's formula. In close 
analogy with the discrete mechanical systems, the curves in the {u/uiy vs. 
N/Ni plane are represented by straight lines. Also in this case, the coordinates 
of the points along these fines provide the critical conditions leading to the 



15 




Figure 9. Nondimensional frequencies vs. nondimensional axial force for the contin- 
uous system in Fig.8. 

system instability in terms of frequency of the excitation and magnitude of 
the apphed compressive axial force. 



3.2 Oscillations and lateral-torsional buckling of beams 



Let us consider a beam of thin rectangular cross-section, constrained at the 
ends so that rotation about the longitudinal axis Z is prevented. Let this beam 
be subjected to uniform bending by means of the application at the ends of 
two moments m contained in the plane YZ of greater flexural rigidity (see 
Fig.lO). 

Considering a deformed configuration of the beam, with defiection thereof in 
the XZ plane of smaller flexural rigidity, and simultaneous torsion about the 
Z axis (sec Fig.lO), bending-torsional out-of- plane vibrations of the beam are 
described by the following partial differential equations: 

where u(z,t) and Lpz{z,t) are, respectively, the out-of-planc deflection and 
the twist angle of the beam cross-section; Ely and Git are the bending and 

torsional rigidities; 11 is the mass of the beam per unit length, and p = Ip/A 
is the polar radius of inertia of the beam cross-section. 
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(d) 



Figure 10. Scheme of a beam subjected to lateral-torsional buckling. 

A solution to the system (47) can be found in the following variable-separable 
form [25]: 



u{z,t) =U{t)r]{z), 
^z{z,t) =<^>{t)iP{z), 



(48) 



where the functions ri{z) and il^{z) are such that the boundary conditions 
rj{0) = r]{l) — T] {0) — r] (l) = ip{0) = — are satisfied. According to 



n'KZ 



Bolotin [25], we can assume ri{z) = ip{z) = sin — — , with n being a natural 

t 

number. In this case, we obtain the following matrix form: 
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GL 
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]l 
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(49) 
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which can be symbohcally rewritten as: 



[M] {q} + [K] {q} - m [Kg] {q} = {0}, (50) 

where {q} = {U, $)^. The mass matrix, [M], the elastic stiffness matrix [K], 
and the geometric stiffness matrix [Kg] in Eq.(50) can be defined in comparison 
with Eq.(49). Looking for a general solution in the form {q} = {qo}e^'^^, we 
obtain: 

{[K] - m [Kg] - [M]) {go} = {0}. (51) 

A nontrivial solution to Eq.(51) exists if and only if the determinant of the re- 
sultant coefficient matrix of the vector {go} vanishes. This yields the following 
generalized eigenvalue problem: 

det (^[K] - m [Kg] - a;^ [M]) = 0, (52) 
where m and cu^ are the eigenvalues of the system. 

As limit cases, if = 0, then we obtain the critical bending moments given 
by the Prandtl's formula: 

m„e = —^ElyGIt, (53) 

whereas, if m = 0, then we obtain the natural flexural and torsional eigenfre- 
quencies of the beam: 



n 



n7r\2 I EL 



, flex / ^ 



I ) \ \^ ' (54) 



Considering a rectangular beam with a depth to span ratio of 1/3 and with 
a thickness to depth ratio of 1/10, the evolution of the first two flexural and 
torsional eigenfrequencies of the system are shown in Fig. 11 as functions of the 
applied bending moment. In this case, the curves in the nondimensional plane 
[LojiJ^^^'^ VS. m/iTLic are no longer straight lines. This fact can be ascribed to 
the coupling between torsional and flexural vibrations of the beam. Moreover, 
when m is increased from zero (pure resonance instability) up to the critical 
bending moment computed according to the Parandtl's formula (pure buckling 
instability), mic, we note that the resonance frequencies related to flexural 
oscillations progressively decrease from u^'^^ down to zero in correspondence 
of the critical bending moments given by the Prandtl's formula. Conversely, 
the resonance frequencies related to torsional oscillations increase. From the 
mathematical point of view, this is the result of the fact that the sum of the 
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Figure 11. Nondimensional flexural and torsional frequencies vs. nondimensional 
bending moments for the continuous system in Fig. 10. 

squares of the two eigenfrequencies for a given value of n is constant when the 
applied moment m is varied. 



4 Finite elements procedure 



When the mechanical system cannot be reduced to the schemes previously 
analyzed, it is possible to apply the Finite Element Method [26,27]. According 
to this approach, the equations of motion for an elastic system with a finite 
number of degrees of freedom can be expressed in matrix form, taking also 
into account the effect of the geometric nonlinearity through the geometric 
stiffness matrix [28,29]. 

For the sake of generahty, let the elastic domain V be divided into subdomains 
Ve, and let each element contain m nodal points, each one having g degrees 
of freedom. In compact form, the displacement vector field defined on the 
element Vg may be represented as: 

{Ve} = [^e] {Se} , 

(gxl) gx{gxm) {gxm)xl ^ ' 

where [r/e] is the matrix collecting the shape functions and {5e} is the nodal 
displacements vector. 
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The deformation characteristic vector is obtained by derivation: 



{9e} = [d] {Ve} , 
(dxl) {dxg){gxl) 



(56) 



where [d] is the kinematic operator relating strains to displacements, whereas 
d denotes the dimension of the strain characteristic vector, e.g. d = 3 for a 
beam in plane or d = 6 for a beam in space. Hence, introducing Eq.(55) into 
Eq.(56), we obtain: 



{Qe} = [d] [Ve] {Se} = [Be] {Se} , 

(dxl) {dxg) gx{gxm) {gxm)xl dx{gxm) {gxm)xl 

where the matrix [Bg] is calculated by derivation of the shape functions. 



(57) 



According to these definitions, we can obtain the following matrix equation 
for each finite element (see [1], Ch.ll for more details): 



[Me] {5e} + {[Ke] - [K,e]) {Se} ^ {0} , 



(58) 



where [Me], [Kg] and [Kg^] denote, respectively, the local mass matrix, the 
local elastic stiffness matrix and the local geometric stiffness matrix of the 
finite element. 

As usual, the local mass and elastic stiffness matrices are given by [26,27]: 
M = / [Ve]''[fA[Ve]dV, (59a) 
[Ke] = I [BeY [H] [Be] dV. (59b) 

The local geometric stiffness matrix can be computed as follows [28]: 

[Kge]^ I [Ge]^ [Se][Ge]dV, (60) 

where the matrix [^'e] is related to the components of the stress field: 



[^e 



[h] Txy [h] Txz [h] 
Txy[h] (yy[h] ryz[h] 
Txz [h] V [h] CTz [h] 



(61) 



where [I^ denotes a unit matrix with dimensions (3x3). The matrix [Gg] is 
related to the first derivative of the shape functions through the differential 



operator 



d 



[Ge] = [d\ [Ve] 

{gxg)x{gxm) {gxg)xg gxigxm) 



(62) 
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where 



d 



- d_ 

f 

- dz 



[h] 



[h] 



(63) 



According to this formulation, we note that the geometric stiffness matrix 
is a function of the stress components through the matrix \S^. In the case 
of a compressive stress field, the geometric stiffness terms become negative 
and reduce the corresponding elements of the local elastic stiffness matrix, 
just as shown for the discrete mechanical systems. We also remark that this 
formulation is quite general, since the information related to the finite ele- 
ment topology is simply included in the matrix [r/e] which collects the shape 
functions and in the differential operator [d] (see [1] for more details). 

By performing the usual operations of rotation, expansion and assemblage 
of the mass, elastic stiffness and geometric stiffness matrices of the element, 
Eq.(58) can be written in global form: 



Looking for the solution to Eq.(64) in the general form {5} = {^oje^'^*, where 
uj is the natural frequency of the system, we can formulate the generalized 
eigenproblem as in the cases discussed above: 



Therefore, the numerical procedure for the determination of the frequency- 
loading multiplier diagram consists in the following steps. 

(1) For a given loading configuration defined by the loading multiplier A, 
determine the stress field according to a linear elastic stress analysis. 

(2) Compute the local mass matrix, the local elastic stiffness matrix and the 
local geometric stiffness matrix for each finite element. 

(3) Perform the rotation, expansion and assemblage operations to obtain the 
global matrices. 

(4) Solve the generalized eigenvalue problem of Eq.(65) and find the eigen- 
frequencies of the system, a;^, with i — 1, . . . , g x n. 

(5) Iterate the above-described procedure for different values of A. 




(64) 




(65) 



21 



5 Discussion and conclusions 



The problems of elastic instability (buckling) and dynamic instability (reso- 
nance) have been the subject of extensive investigation and have received a 
large attention from the structural mechanics community. Nonetheless, the 
study of the interaction between these elementary forms of instability is still 
an open point. 

The phenomenon of flutter instability of the Tacoma Narrows Bridge occurred 
on November 7, 1940, can be reinterpreted as the result of such a catastrophic 
interaction. This cable-suspended bridge was solidly built, with girders of car- 
bon steel anchored in huge blocks of concrete and was the first of its type to 
employ plate girders to support the roadbed. While in the earlier designs any 
wind would simply pass through the truss, in the new design of the 1940 the 
wind would be diverted above and below the structure. Shortly after construc- 
tion, it was discovered that the bridge would sway and buckle dangerously in 
windy conditions. This resonance was flexural, meaning the bridge buckled 
along its length, with the roadbed alternately raised and depressed in certain 
locations. However, the failure of the bridge occurred when a never-before-seen 
twisting mode occurred. 

A Report to the Federal Works Agency [30] excluded the phenomenon of 
pure forced resonance as the actual reason of instability: " ...it is very improb- 
able that resonance with alternating vortices plays an important role in the 
oscillations of suspension bridges. First, it was found that there is no sharp 
correlation between wind velocity and oscillation frequency such as is required 
in case of resonance with vortices whose frequency depends on the wind ve- 
locity...'" . A new theory for the interpretation of these complex aerodynamic 
instabilities was developed by Scanlan [21, 22] and then elaborated by var- 
ious researchers [31-34]. Basically, the so-called flutter theory considers the 
following equation of motion for the mechanical system in the finite element 
framework [9,32]: 

[M] {5} + [C] {6} + [K] {5} = {F}„, + {F}^, , (66) 

where and {F}^^ are, respectively, the motion-independent wind force 

vector and the motion-dependent aeroelastic force vector. A special attention 
is given to the structural damping, which is included in the equations of motion 
through the damping matrix [C]. The motion-dependent force vector is then 
put in relationship with the nodal displacements of the system, {5}, and the 
nodal velocities, {d}, according to the flutter derivative matrices, [K*] and 
[C*], that are empirically determined in the wind tunnel by using section 
models of the bridge. As a result, the problem becomes highly nonlinear, and 
the flutter velocity, Ucr, and the flutter frequency, cjcr, can be determined from 
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the following eigenproblem: 

det (^[K] - ^pUl [K*] - [M] + [C] - ^pt/e.^cr [C*]) = 0, (67) 

where p is the air density and IjlpUl^ is the wind pressure. 

It is important to remark that this eigenproblem shares most of the fea- 
tures of the generalized eigenproblem that we have analyzed in the present 
study. In fact, in both cases, two eigenvalues have to be determined from 
the eigenequation. However, the flutter theory gives prominence to the role 
played by the structural damping, although being generally less than 1% (see 
e.g. [32]). Moreover, the value of the mechanical damping seems to represent 
a sort of free parameter in the model. In fact, this parameter is usually as- 
sumed, like in [32], rather than experimentally evaluated. Another difference 
with our proposed approach relies in the geometric stiffness matrix, which 
is not taken into account in the current flutter theory. However, due to the 
large structural displacements, the flutter derivative matrix \K*\ plays a very 
similar role, although being experimentally obtained, rather than analytically 
computed. 

In conclusion, it seems to be possible to reinterpret the phenomenon of aeroe- 
lastic instability as the result of the interaction between pure resonance and 
pure buckling instabilities. According to our approach, the geometric stiffness 
matrix has a preeminent role and the mechanical damping can be neglected, as 
usually done in most of the structural engineering applications. On this line, 
the collapse of the Tacoma Narrows bridge can be considered as the result of 
the interaction between buckling (related to the wind pressure proportional 
to the square of the wind velocity) and resonance (caused by the frequency of 
the wind gusts). Thus, this would give a new explanation on why the Tacoma 
Narrows bridge failure took place under moderate wind velocities (wind pres- 
sure lower than the critical buckling load) and wind gusts frequencies different 
from the natural frequencies of the bridge. Future developments of the present 
work will regard the assessment of the proposed approach to the analysis of 
bridge instabilities, as well as the comparison with the classical flutter theory 
on the basis of real case histories. 



6 Appendix: orthonormality of the eigenfunctions of deflected beams 
subjected to an axial force 

As is well-known, the eigenfunctions rii of deflected beams computed according 
to pure modal analysis are orthonormal functions. It is possible to demonstrate 
that this fundamental property still holds when the beam is subjected to an 
axial load, A^, as that shown in Fig.8. We may in fact write Eq.(38b) for two 
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different eigensolutions: 



alrik. 



(68a) 
(68b) 



Multiplying Eq.(68a) by rjk and Eq.(68b) by rjj, and integrating over the beam 
length, we obtain: 



VjVk^dz + 0^ VjVk^z ^Oilj^ VjVkdz. 



(69a) 
(69b) 



Integrating by parts the left-hand sides, the foregoing equations transform as 
follows: 



(70a) 





I 








Vkvr[ 










/ 

Jo 



rl rl 
Jo ^^^^^^^ ^ Jo ^kVA^, 

1 + J\'/Vkdz + P'[7j,7j,i 



vWk 



rl rl 



(70b) 



When each of the two ends of the beam is constrained by a built-in support 
[f] = Tj^ = 0), or by a hinge {rj = rj^^ = 0), the quantities in square brackets 
vanish. On the other hand, when the end in ^ = is either unconstrained 
{rj^^^ — rj^^ = 0), or constrained by a double rod {rj^^^ — rj^ — 0), the remaining 
end of the beam has to be constrained either by a built-in support (r) = r)^ ~ 
0), or by a simple support [r] = rj^^ = 0). For both configurations, only the 
terms [r]ir]k\Q are different from zero. 

In any case, subtracting member by member, these quantities are canceled 
and we have: 



(a| - al) rij7]k = 0, 
which leads to the orthonormality condition: 

/ VjVk = Sij, 

J 



(71) 



(72) 



where 5ij is the Kronecker delta. Thus, when the eigenvalues are distinct, the 
integral of the product of the corresponding eigenfunctions vanishes. When, 
instead, the indices j and k coincide, the condition of normality reminds us 
that the eigenfunctions are defined neglecting a factor of proportionality. 
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